function matrix = sparse_pentadiagonal(main,p1,m1,p2,m2,m)
% construct sparse pentadiagonal matrix with zero boundary
% condition compensated on the first off-diagonals

a = zeros(m^2,m^2);

p1v = repmat(p1,m^2,1);
m1v = repmat(m1,m^2,1);
mainv = repmat(main,m^2,1);
p2v = repmat(p2,m^2,1);
m2v = repmat(m2,m^2,1);

% columns of b represent diagonal bands in constructed matrix
b = [m2v m1v mainv p1v p2v]; 
d = [-m, -1, 0, 1, m];       % diagonal band location specifier
matrix = spdiags(b,d,a);

% insert zeros in first off-diagonals to compensate for boundary condition

for i=2:m^2-1
  mod4 = mod(i,4);
  if mod4 == 0
    matrix(i,i+1)=0;
  elseif mod4 == 1
    matrix(i,i-1)=0;
  end
end
